Create a table summarizing data layers within each MTBS fire. Layers were assembled and reduced in Google Earth Engine, then the summary file exported to Google Drive. Downloaded the summary table locally from Google Drive, then ran this script to visualize and explore the data.

# how many fires in this data set? 
message("# of fires in dataset: ")
## # of fires in dataset:
message(nrow(fire_df))
## 37
# create output directory for figures
if(!dir.exists(here::here("figures"))){
  dir.create(here::here("figures"))
}
# create a table to summarize attributes of interest per fire 
fire_table <- as.data.frame(fire_df) %>% 
  # select columns of interest
  dplyr::select(Fire_Name, Year, Acres, gedi_coverage, lodgepole, ponderosa, spruceFir, disturbed_burned, disturbed_unspecific, 
                #disturbed_logged, regenerating_disturbed, regenerating_harvested
                ) %>%
  # reorder the rows based on multiple variables 
  dplyr::arrange(desc(gedi_coverage), desc(lodgepole), desc(ponderosa), desc(spruceFir), desc(disturbed_burned)) %>%
  # rename the columns for interpretation
  dplyr::rename("Fire Name" = Fire_Name, 
                "Year" = Year,
                "Acres" = Acres,
                "GEDI %" = gedi_coverage,
                "Lodgepole %" = lodgepole,
                "Ponderosa %" = ponderosa,
                "Spruce Fir %" = spruceFir,
                "Disturbed Burned %" = disturbed_burned,
                "Disturbed Unspecified %" = disturbed_unspecific,
                #"Disturbed Logged %" = disturbed_logged,
                #"Regenerating Disturbed %" = regenerating_disturbed,
                #"Regenerating Harvested %" = regenerating_harvested
                ) 

# add a categorical column with dominant forest type 
temp <- fire_df %>% 
  dplyr::select(Fire_Name, lodgepole, ponderosa, spruceFir) %>% 
  tidyr::gather(major_forest_type, forest_percent, lodgepole:spruceFir) %>% 
  dplyr::group_by(Fire_Name) %>% 
  dplyr::slice(which.max(forest_percent)) %>% 
  dplyr::select(Fire_Name, major_forest_type) %>% 
  st_set_geometry(NULL)
fire_df <- fire_df %>% 
  dplyr::left_join(temp)
remove(temp)



kableExtra::kable(fire_table) %>%
  kableExtra::kable_styling(bootstrap_options = "striped", 
                            full_width = F, 
                            position = "left")
Fire Name Year Acres GEDI % Lodgepole % Ponderosa % Spruce Fir % Disturbed Burned % Disturbed Unspecified %
MAD CREEK 2001 1435 99.7 1.1 0.0 74.6 0 0
BIG ELK 2002 4347 98.5 19.3 47.9 0.0 0 0
CRYSTALFIRE 2011 2855 95.3 2.7 54.2 0.0 0 0
LOWER NORTH FORK FIRE 2012 3436 93.8 0.2 31.2 0.0 0 0
#6 1989 1547 93.2 0.0 6.7 0.0 0 0
SNAKING 2002 2080 91.4 0.2 65.4 2.7 0 0
FOURMILE CANYON 2010 5865 91.1 1.2 67.0 0.0 0 0
BIG RED 2017 2868 90.0 80.6 0.0 0.9 0 0
UTE CREEK 1994 3225 89.2 2.1 0.0 36.0 0 0
PICNIC ROCK 2004 9014 88.2 0.0 15.9 0.0 0 0
MT. ZIRKEL COMPLEX (HINMAN) 2002 17627 86.4 10.1 0.1 23.1 0 0
BUFFALO CREEK 1996 9792 85.5 0.0 9.4 0.0 0 0
OLD MAN 2001 1169 83.4 0.0 0.0 0.0 0 0
MEADOW CREEK 2010 1432 83.3 3.8 0.1 8.9 0 0
MT. ZIRKEL COMPLEX (BURN RIDGE) 2002 15478 71.8 4.4 0.0 58.0 0 0
BOBCAT 2000 9131 69.8 1.4 75.9 0.0 0 0
BIG FISH 2002 17273 68.3 3.2 0.0 58.2 0 0
HAYDEN PASS 2016 17978 61.7 2.0 21.5 15.0 0 0
BEAVER CREEK 2016 44222 60.4 47.7 0.0 3.6 0 0
WALDO CANYON 2012 20112 58.6 0.3 47.2 0.2 0 0
LOST SOLAR 2016 5142 56.3 3.1 0.0 52.9 0 0
SUNNYSIDE 1989 1416 47.8 0.0 27.5 0.0 0 0
CANYON 1988 1092 47.5 0.0 40.8 0.0 0 0
HIGH PARK 2012 90769 46.8 13.9 49.6 1.1 0 0
LOST LAKES 2002 5888 44.9 0.2 0.0 51.8 0 0
UNNAMED 1989 1305 43.6 0.0 43.7 0.0 0 0
DUCKETT 2011 4705 41.2 0.4 17.4 0.0 0 0
COAL SEAM 2002 12004 37.3 0.1 0.0 0.5 0 0
OVERLAND 2003 3232 35.6 1.8 79.6 0.0 0 0
HAYMAN 2002 129417 25.4 0.9 57.4 0.6 0 0
GREEN CREEK 2002 4342 21.6 16.4 0.0 55.5 0 0
SPRINGER 2012 1666 14.7 0.0 83.9 0.2 0 0
RED FEATHER 2017 1582 13.0 0.4 79.5 0.0 0 0
GRACE CREEK 1988 1455 11.7 32.1 0.5 3.4 0 0
SPRING CREEK 2002 11668 6.5 0.7 0.0 45.5 0 0
HIGH MEADOWS 2000 9607 0.0 0.0 75.0 0.0 0 0
SCHOONOVER 2002 2821 0.0 0.0 73.2 0.0 0 0
# count number of fires with each dominant forest type 
kableExtra::kable(table(fire_df$major_forest_type),
                  col.names = c("Dominant forest type", "# of fire events")) %>%
  kableExtra::kable_styling(bootstrap_options = "striped", 
                            full_width = F, 
                            position = "left")
Dominant forest type # of fire events
lodgepole 4
ponderosa 22
spruceFir 11

Explore the fire data

Map

library(leaflet)
library(rgdal)
library(geojsonio)

# read geojson file for mapping with leaflet
fire_gjson <- rgdal::readOGR(fire_filename)
## OGR data source with driver: GeoJSON 
## Source: "C:\Users\tyler\Documents\Field Planning\CAREER\data\fires_stats_20210414.geojson", layer: "fires_stats_20210414"
## with 37 features
## It has 14 fields
# add a column with fire perimeter centroid point locations for map labels
fire_df$centroid_lon <- NA
fire_df$centroid_lat <- NA
fire_df$centroid_label <- NA
for (i in 1:nrow(fire_df)){
  # create popup label using fire name and year
  content <- paste(sep = "<br/>",
    fire_df$Fire_Name[i], fire_df$Year[i])
  #print(content)
  fire_df$centroid_label[i] <- content
  
  # get lon, lat coordinates of geometry centroid  
  fire_df$centroid_lon[i] <- st_centroid(fire_df$geometry[i])[[1]][1]
  fire_df$centroid_lat[i] <- st_centroid(fire_df$geometry[i])[[1]][2]
}


# color palette
pal <- colorNumeric(
  palette = "viridis",
  domain = fire_gjson$Year
)

# map the fire polygons
map <- leaflet(fire_gjson) %>%
  addTiles() %>%
  addPolygons(stroke = FALSE, smoothFactor = 0.2, fillOpacity = 0.7,
    color = ~pal(Year), group = "MTBS fires"
  ) %>%
  addLegend("bottomright", pal = pal, values = ~Year,
    title = "Fire Year",
    # get rid of the comma in years (1,990 becomes 1990) for legend entries
    labFormat = labelFormat(big.mark = ''),
    opacity = 1
  ) %>% 
  # set base map type
  addProviderTiles(providers$CartoDB.Positron, group = "Carto") %>%
  addProviderTiles(providers$Stamen.Terrain, group = "Stamen")

# Icons 
icons <- makeAwesomeIcon(
  icon = "fire",
  iconColor = "gray",
  library = "fa",
  markerColor = "gray"
)


# Add Fire Name popup labels to each polygon centroid
  map <- map %>%
    #addMarkers(lng = fire_df$centroid_lon,
    addAwesomeMarkers(lng = fire_df$centroid_lon,
               lat = fire_df$centroid_lat,
               icon=icons,
               popup = fire_df$centroid_label) %>% 
    addLayersControl(baseGroups = c("Carto", "Stamen"), 
                   overlayGroups = c("MTBS fires"))  


# display map
map

Fire timeline

# calculate minimum and maximum fire years
year_min <- min(fire_df$Year)
year_max <- max(fire_df$Year)

# set forest type colors for the figures
                       # lodgepole  ponderosa  spruceFir
#forest_type_colors <- c("#005a32", "#74c476", "#e5f5e0") 
forest_type_colors <- c("#04663B", "#41ab5d", "#c7e9c0")


# reorder the rows based on year
fire_df <- fire_df %>% 
  dplyr::arrange(desc(Year)) 

# calculate point size - scaled relatively by fire size
tmp <- log(fire_df$Acres / max(fire_df$Acres))
tmp <- tmp + abs(min(tmp)) + 1.25

# create timeline figure 
ggplot(fire_df) + 
  geom_point(aes(x = reorder(Fire_Name, Year), y = Year, 
                 # color points based on dominant forest type
                 fill = major_forest_type), 
             # add a black outline around each point
             color = "black", shape=21, 
             # size of each point, thickness of outline
             size = tmp, 
             stroke = 0.3) + 
  scale_y_continuous(breaks = seq(year_min, year_max, by = 4)) + 
  labs(title = "Fire Event Timeline", y = "Year", x = "Fire Name\n",
       fill = "Dominant forest type") + 
  # Put fire name on Y axis, year on X axis
  coord_flip() + 
  # set the point fill colors using hex codes 
  scale_fill_manual(values = forest_type_colors,
                    labels = c("Lodgepole pine", "Ponderosa pine", "Spruce/Fir")) + 
  # set the point size in legend 
  guides(fill = guide_legend(override.aes = list(size=4))) + 
  # clean theme for plot with white background
  theme_bw() 

# Save figure to file
ggsave(filename = (here::here(file.path("figures",paste0(out_label,"-fire_timeline.png")))),
       width = 8, height = 6)

Fire year histogram

hist_year_title <- paste("Histogram: Fire Years,", year_min, "-", year_max)

# create histogram, fire count per year.
# color the bars based on the dominant forest type of each fire event. 
ggplot(fire_df, aes(x=Year, fill = major_forest_type)) + 
  geom_histogram(binwidth=1, color="black", 
                 # set the outline thickness and bar transparency. 
                 size = 0.2, alpha=0.9) + 
  labs(title = hist_year_title, y = "Count", fill = "Dominant forest type") + 
  # x axis label years from min to max year in increments of 4
  scale_x_continuous(breaks = seq(year_min, year_max, by = 4)) + 
  # set color scale of the dominant forest types
  scale_fill_manual(values= forest_type_colors,
                    labels = c("Lodgepole pine", "Ponderosa pine", "Spruce/Fir")) + 
  theme_bw()

# Save figure to file
ggsave(filename = (here::here(file.path("figures",paste0(out_label,"-fire_hist_years.png")))),
       width = 7, height = 5)

Fire size histogram

# color the bars based on the dominant forest type of each fire event. 
           # convert from Acres to Hectares? * 0.404686
ggplot(fire_df, aes(x = Acres * 0.404686)) + 
  geom_histogram(aes(fill = major_forest_type), color="black", 
                 # set the outline thickness and bar transparency. 
                 size = 0.2, alpha=0.9) + 
  labs(title = "Histogram: Fire Size", y = "Count", x = "Size [Hectares]",
       fill = "Dominant forest type") + 
  # set color scale of the dominant forest types
  scale_fill_manual(values = forest_type_colors,
                    labels = c("Lodgepole pine", "Ponderosa pine", "Spruce/Fir")) + 
  theme_bw()

# Save figure to file
ggsave(filename = (here::here(file.path("figures",paste0(out_label,"-fire_hist_size.png")))),
       width = 7, height = 5)

Fire size as a function of year

ggplot(fire_df, aes(x=Year, y=Acres * 0.404686 , 
                    log="y", label = Fire_Name)) + 
  geom_point(aes(# color points based on dominant forest type
                 fill = major_forest_type), 
             # add a black outline around each point
             color = "black", shape=21, 
             # size of each point, thickness of outline
             size = 4, stroke = 0.5) + 
  labs(title = "Fire size vs year", 
       y = "Hectares [log transformed]", 
       fill = "Dominant forest type") + 
  # log transform the y axis to space out the points more 
  scale_y_continuous(trans='log2') + 
  scale_x_continuous(breaks = seq(year_min, year_max, by = 4)) + 
  scale_fill_manual(values = forest_type_colors,
                    labels = c("Lodgepole pine", "Ponderosa pine", "Spruce/Fir")) + 
  # add Fire Name labels
  geom_label_repel(aes(label = Fire_Name),
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   segment.color = 'grey50') + 
  
  theme_bw() + 
  # font size
  theme(axis.text.y = element_text(size=16),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=16),
        axis.title.x = element_text(size=16),
        # move the legend to the top of the figure
        legend.position="top") 

# Save figure to file
ggsave(filename = (here::here(file.path("figures",paste0(out_label,"-fire_size_vs_year_scatterplot.png")))),
       width = 7, height = 9)